Topological insulators in strained graphene at weak interaction 
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The nature of the electronic ground states in strained undoped graphene at weak interaction 
between electrons is discussed. After providing a lattice realization of the strain-induced axial mag- 
netic field we numerically find the self-consistent solution for the time reversal symmetry breaking 
quantum anomalous Hall order-parameter, at weak second-nearest-neighbor repulsion between spin- 
less fermions. The anomalous Hall state is obtained in both uniform and nonuniform axial magnetic 
fields, with the spatial profile of the order-parameter resembling that of the axial field itself. When 
the electron spin is included, the time reversal symmetric anomalous spin Hall state becomes slightly 
preferred energetically at half filling, but the additional anomalous Hall component should develop 
at a finite doping. 

PACS numbers: 71.10.Pm, 71.70.Di, 73.22.Pr 
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Graphene, being a flexible membrane, is always wrin- 
kled to certain degree. One can think of the effect of the 
wrinkles on the quasirelativistic Dirac low-energy excita- 
tions in graphene as the random, static component of the 
non-Abelian (or "axial" ) gauge field, which arises from 
modified hoping amplitudes in the tight-binding model, 
and which preserves the time-reversal invariance. Wrin- 
kles are typically randomly distributed, yielding roughly 
zero total flux of the axial field. If one deliberately bulges 
graphene, however, the total axial magnetic flux may be- 
come finite. In a recent experiment 1 graphene was de- 
posited over a metallic substrate and cooled down. The 
mismatch of graphene's and substrate's compressibilities 
subjects the graphene layer then to a strain, resulting 
in a surprisingly uniform axial magnetic field as high as 
~ 350 T. We will argue here that these may be the ideal 
conditions for a possible observation of the dynamically 
induced quantum anomalous Hall (AH) and the quantum 
anomalous spin Hall (SH) states, the former predicted be- 
fore in Ref. |2|. In so doing we will go beyond the usual 
Dirac approximation and present numerical evidence for 
the effect, even in the more general case of a spatially 
non-uniform axial field. 

The gist of the effect lies in the index theorem 3 -, which 
equally well applies to the axial and real magnetic fields: 
a finite arbitrary field, in the continuum, produces single- 
particle states at zero-energy, proportional in number 
to the total flux. This zero-energy degenerate band in 
graphene is half-filled, and it is energetically advanta- 
geous to split it and populate only the lower half. We 
show here that this is what happens in presence of even 
a weak second-nearest-neighbor repulsion, which, when 
strong, is thought to favor the formation of the AH state, 
even without any axial fluxi* The time reversal symme- 
try (TRS) breaking AH order parameter (OP)£ is found 
in our solution to be spatially distributed similarly to 
the axial flux itself, and would thus be close to uniform 
in a uniform field in the experiment. We first present 
the detailed numerical results of our self-consistent cal- 



culation for the spinless fermions. With the restoration 
of spin, other ordered states, such as the spin polarized 
ferromagnetic stated and the anomalous SH insulator 7 
also become possible. Neglecting possible effects of the 
Hubbard on-site interaction we find that, although de- 
generate at the mean-field level, fluctuations in this case 
favor slightly the SH state in the finite axial field and at 
weak coupling. 




FIG. 1: (Color online) x(^) attached to each site; its value 
increases in the order black, red, blue, magenta. Thick bonds 
represent modified hopping amplitude. Sections (a,c, e), and 
(b,d,f) are connected by 2n/3 rotations. 

To construct a lattice implementation of the axial 
field, we recall first the low-energy Dirac Hamiltonian 
in graphene. In the presence of an axial field, it may be 
written as^ 

H[a] = i loli {-idi - i^a^x)) = e x ^°H[0]e^°, 

(1) 

where i = 1,2. Neglecting spin, the Dirac Hamiltonian 
acts on the four component fermion, defined as 

V T (q) = (u(K+q),v(K+q),u(-K+q),v(-K+q)). (2) 

u (v) is the annihilation operator on sublattice A (B), 
and cii(x) = eijdjx{%)- The four-component Hermi- 
tian 7-matrices are 70 = <jq <g> 03 , 71 = (73 01, 72 = 
Co ® 01,73 = o\® 0-2,75 = 02 <g> 02 9 ' 10 . In the pres- 
ence of an axial field, the emergent chiral SU C (2) sym- 
metry of H[0], generated by {73, 75, 27375}, is reduced to 
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U c (l), generated by ij 3 ^gll. The pseudo magnetic field 
is 6(f) = eijd^jix) = d 2 x(x). 

From Eq. ([l} one can write the zero energy states in 
the presence of axial field as 

* ,„ [a] (x) oc e-*<3Tb ^ [0] (3) 

The matrix 70 in the exponent alternates in sign on A and 
B sublattices. Since x(f) is an increasing function of the 
distance, 3 only the sublattice A supports normalizable 
zero-energy states. Zero-energy states on the sublattice 
B will diverge at infinity, i. e., in a finite system would be 
localized near the boundary. This suggest an introduc- 
tion of the axial gauge potential on lattice as the following 
modification of the nearest-neighbor hopping integrals: 

t a p = e x{a) t e~ x{p \ (4) 

where a € A, /3 € B, and t(= 1) is the uniform 
hoppingJ^ 

Upon defining a quantity 7Z, counting the minimal 
number of bonds required to reach a particular site from 
the central hexagon, we assign x(^) to each site, de- 
pending on whether it belongs to the sublattice A or B, 
in the following way: when TZ is odd, x{A) > x(-3)> 
and when it is even, x(A) = x{B). This is presented 
in Fig. [TJ TZ here plays the role of the radial coordi- 
nate, and for all six sites in the central hexagon in Fig. Q] 
TZ = 0, for example. Then along each bond with modi- 
fied hopping amplitude x{A) > x(-S): in agreement with 
Eq. (01 . Such modification leaves the honeycomb lattice 
invariant under the C3 symmetry, and an axial vector 
potential a = (or, a$) ~ (0, dx(TZ) / dTZ) is introduced in 
the system. If x(7£) oc 1Z 2 , the axial flux enclosed by the 
system is roughly proportional to its area; corresponding 
to an approximately uniform axial magnetic field. On 
the other hand, a bell-shaped axial field, localized around 
the center of the system, can be obtained by choosing 
x 1ip:7C,:. 

Let us first take a spinless fermion hopping on a hon- 
eycomb lattice of 2400 sites-i 3 -, subject to uniform and 
nonuniform axial fields. Since the non-interacting Hamil- 
tonian is bipartite, the diagonalization always gives a 
particle-hole symmetric spectrum, with a finite number 
of states localized close to the zero energy. With the 
bond configuration as in the above, we indeed find that 
the near-zero-modes that are on the sublattice A are lo- 
calized in the bulk of the system, while those that are on 
the sublattice B are localized near the system's bound- 
ary. Furthermore, when the field is roughly uniform, iso- 
lated windows of energy W ~ t/10 where the number of 
states increases with the axial flux appear symmetrically 
around zero. These states appear to form the first Lan- 
dau level (LL), and we define the mean energy of these 
states, which scales roughly as \/b, as the first LL energy 
(£1 ) . The lack of perfectly sharp LL quantization in our 
calculation is due to the finite size effects, and also to the 
local modification of the Fermi velocity^. 

Next we turn on the second-nearest-neighbor interac- 
tion (V) among the fermions, still kept spinless. The 
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FIG. 2: (Color online) AH OP (in units of i) on A(left), 
B (right) sublattice in the presence of uniform axial field 
b = 0.0256 o . 6 = h/(ea 2 ) « 10 4 T, is the field associated 
with graphene's lattice spacing a w 2.5A. h is Plank's con- 
stant, e is electronic charge. Interaction strength reads as V = 
1.5,1.27,1.0,0.75,0.5 from top to bottom, r = n(l,2, •••) 
corresponds to nth ring around the center of the system 1 ^. 

interacting Hamiltonian reads 

H = ^2 t a p (c^cp + C^Caj + V ^ n a n P- (5) 

(«,/8) «a,/3)> 

((. . . )) stands for the sum over the nearest-neighbor and 
second-neighbor sites, and n a is the fermion occupation 
number on the site a. The usual Fock decomposition 
yields an effective single-particle Hamiltonian 

HsP = t a p (c^Cp + C^C a j + ^ ^lapCpCa + H.c). 

We will assume the intra-sublattice circulating current 
i]af3 — V(cL c i3) ^° be purely imaginary, and oriented in 
opposite directions on the two sublattices. It will play 
the role of the TRS breaking AH order parameter—, and 
we proceed to determine it self-consistently on a finite 
honeycomb lattice. In the absence of an axial field (x = 
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FIG. 3: (Color online) Left: Difference of the AH OP (in 
units of t) in bulk (rojj) and edge (m^) on A-sublattice as 
a function of the system size for V=1.27(red), l.O(black), 
0.75(blue) with b = 0.015&0- Right: Scaling of anoma- 
lous Hall OP with axial magnetic fields (b/bo). The red, 
black, magenta and green dots respectively correspond to 
V = 1.5, 1.27(Vb), 1.0, 0.75.0.5. 

0) a non-zero self-consistent solution of the AH order 77 
on a honeycomb lattice of 600 sites is found only for V > 
1.27—. The amplitudes of the OP on both sublattices 
are then equal. For V < V c , the self-consistent value of 77 
vanishes everywhere in the system. Hence, we can define 
V c = 1-27 as the zero-field critical interaction. This value 
is very close to the one found analytically^^ 
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After the introduction of the uniform axial field in the 
same system, we search for the self-consistent solution 
of r\ at sub-critical interactions V < V c as well. Typical 
distributions of the OPs then are shown in Fig. [5J For 

V < V c , the TRS breaking OP on the sublattice A (B) 
clearly forms in the bulk (boundary) of the system. This 
spatial separation of the OPs on two sublattices follows 
the structure of the near-zero-energy states. The OP 
on the sublattice A is roughly uniform in the bulk (for 
r < 6), and with the inhomogeneity disappearing with 
increase of the system size (see Fig.[3j left). As the inter- 
action gets stronger (V ~ V c ), the effect of the axial field 
becomes irrelevant and OPs on two sub-lattices become 
comparable; see the red, black, and blue curves in Fig. 
[21 The effective single-particle Hamiltonian Hsp in Eq. 
(O preserves the C3 symmetry of the bond configuration 
in Fig. 1. 

The OP on the sublattice A (averaged up to r = 6) 
scales linearly with the uniform axial field (b) when 

V <C V c (Fig. |31 right). The scaling becomes sub-linear 
for intermediate strength of the interaction. At the zero- 
field criticality (V = 1.27), the gap appears to scale as 
Vo. This behavior can also be confirmed from the lin- 
ear scaling of m 2 with b which passes though the origin, 
as shown in Fig. |4] (left). Finally, for V > V c the OP 
saturates to a finite value at zero field. With our defini- 
tion of the first LL energy E\ , the universal ratio of E\ 
to the mass gap at V = Vc is found here to be ~ 5.78. 
The scaling of the TRS-breaking OP is thus similar to 
the scaling of the chiral-symmetry-breaking OP with the 
real magnetic field found previousl y 17 ! 18 . The universal 
ratio is also reasonably close to the value 5.985, obtained 
analytically for the chiral symmetry breaking OP. 18 

Existence of the near zero energy states even when the 
axial magnetic field assumes a spatially non-uniform pro- 
file, allows the formation of the AH order for V -C V c . 
With a bell-shaped field, localized around the center of 
the system, a typical distribution of the OP on the sub- 
lattice A, computed self-consistently on a 726 site hon- 
eycomb lattice, is shown in the right of Fig. S) For 

V <C V c , the OP forms only in the vicinity of the lo- 
calized flux. This behavior of the OP on A sub-lattice is 
generic and the OP disappears towards the boundary of 
larger systems^. However, for sub-critical interactions 
the OP is also accompanied by a real amplitude, which 
appears to be a finite size effect. On the other hand, 
when V ~ V c , the OP starts to become considerable ev- 
erywhere in the system. Behavior of the OP on the sub- 
lattice B is qualitatively similar to that in the presence 
of uniform axial field. 

Let us now restore the spin degrees of freedom, and 
define the 8-component spinor W = (^-j-, where 
(7 =r, I are the spin projections along the z-axis. The 4- 
componcnt spinors 'f a for each spin projection take the 
form of \& in Eq. ([2]) . The Dirac Hamiltonian in this basis 
is H]j[a] = (T <E) H[a\. The finite range components of the 
Coulomb interaction then allow the formation of various 
other ordered phases besides the AH insulator; for exam- 
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FIG. 4: (Color online) Left: Linear scaling of m 2 with b for 
V = V c = 1.27. Right: TRS breaking OP for V=1.5(green), 
1.27(magenta), l.O(black), 0.75(blue) with total flux 7.85<!>o 
of a bell-shaped non-uniform axial field. $0 = h/e is flux 
quantum. 

pie, the anomalous SH insulator, or the spin-polarized 
ferromagnetic state. For example, the on-site Hubbard 
repulsion U typically favors the spin-polarized stated. On 
the other hand, the second-nearest- neighbor repulsion 
prefers the AH or the SH insulators^i Depending on 
the relative strengths of the finite ranged components 
of the Coulomb interaction there will be a competition 
among various ordered states in strained graphene. In 
this work, however, we will take only the second-nearest- 
neighbor repulsion into account, and study the competi- 
tion between the AH and SH insulators. Possible appear- 
ance of the ferromagnetic ground state for strong on-site 
Hubbard U, and its competition with the AH and the SH 
orders will be addressed in a separate publication. 




FIG. 5: (Color online) Left: Flow of the couplings g c (red), 
and gh (black). Right: Proposed splitting of the zero energy 
subspace with both SH and AH order present. All the states 
are localized on one sub-lattice. Both SH and AH order pa- 
rameter are formed when the chemical potential is at /xi. 

The SH and the AH OPs arc C = (* f [a <8> 17172] *), 
Co = [<Jo £3> 47172] respectively. The SH insulator 
is even under the TRS, which now also includes the usual 
reversal of the spin, whereas the AH state is oddJ^ Notice 
that the matrices appearing in both the AH and SH OPs 
anticommute with the Dirac Hamiltonian Ho[a\. Hence, 
both of these two OPs upon developing a finite expecta- 
tion value would not only split the states at zero energy, 
but would also shift downward the occupied states with 
negative energies. For example, if the axial field is uni- 
form, the LLs at ±\/2nb get pushed to ±V2n& + X 2 , 
where X = Cq or \C\. Hence, at half-filling when all the 
states at negative (positive) energies are filled (empty), 
it is energetically highly advantageous to develop such a 
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mass-gap . 

Let us now address the competition between the AH 
and the SH insulators in the following way. Consider first 
the interaction Lagrangian in the continuum, and at zero 
field, as: 

Lint =gh (*V «)i7i72*) 2 +5 c (* t o ; ®«7i72^) 2 - (7) 

The interactions gh < and g c < favor the AH and SH 
insulators, respectively. In graphene at b = any weak 
electron-electron interaction is irrelevant, and to place 
the system in an ordered phase the interaction needs 
to be sufficiently strong^iS. Assuming only the second- 
nearest- neighbor repulsion V > 0, gh, g c ~ — V at the lat- 
tice scale, and they are equally irrelevant if V/t <C 1. At 
the mean-field level therefore the AH and the SH states 
are degenerate. As one integrates out the fast Fourier 
modes in momentum-shell A/s < \k\ < A, where s > 1, 
however, these two couplings flow to zero differently. This 
is a consequence of the fact that the corresponding OPs 
obviously break different symmetries, one of which dis- 
crete, and the other continuous. To the second order, the 
flow of these two couplings is given by the /3-functions 
{fix =dX/d\ogs) 

f3 9h = -9h{l + igh - 3.9c), /3g c = -g c (l + 5g c - g h ), (8) 

after rescaling {2A/n 2 )gh yC — > gh,c~ 

Negative linear terms in the /3-functions imply that any 
spontaneous symmetry breaking in graphene at zero field 
can occur only at strong interactions. In the presence of 
an axial or standard magnetic field, however, the order- 
ing will occur even at weak interactions, due to the zero 
energy states, as we have argued above. When the ax- 
ial field is roughly uniform we will choose the final value 
of the parameter s ~ lb /a 3> 1, where a is lattice spac- 
ing and lb ~ and determine the effective values of 
the couplings at the scale of the axial magnetic length. 
We therefore numerically solve the coupled flow equa- 
tions with g c = gt, — V when s = 1. The less irrelevant 
coupling at the scale s ~ lb /a ^> I will then give rise to 
the dominant instability in strained graphene. The flows 
of the coupling constants gh and g c are shown in Fig. [5] 



(left), and one can see that g c is less irrelevant than gh 
for any s > 1. The leading instability at weak coupling 
in the presence of the axial field is therefore the SH state. 
The same outcome is also found at strong coupling and 
at zero field, in agreement with the previous study. 4 

When the chemical potential is close to the first ex- 
cited states at ±|G| the system can develop a gap by 
breaking the TRS and by developing an additional AH 
order, as shown in Fig. [5] (right). The resulting single- 
particle excitation gap of the AH state is then 2 Co, and 
the zero-energy LL is 1/4 or 3/4 full. The scaling of 
the AH OP, which one can possibly realize in strained 
graphene in this way and away from the neutrality point 
is qualitatively similar to the one we have computed nu- 
merically at the neutrality point. When the axial field 
is non-uniform the same mechanism still operates, with 
the AH OP then developing mostly in the vicinity of the 
localized flux at weak interaction. 

To summarize, we proposed a specific modulation of 
the nearest-neighbor hopping amplitude that captures 
the effect of time reversal symmetric axial magnetic 
fields. Although such a field can also be spatially nonuni- 
form, it always produces a finite number of states near 
zero energy. We show that such a spectrum is conducive 
to formation of topologically non-trivial ordered phases 
such as the AH and SH insulators, even for weak repulsive 
interactions. We provide numerical evidence that spinless 
fermions can be in the AH phase even at weak second- 
nearest-neighbor repulsive interaction, with the magni- 
tude of the OP depending linearly on the axial field. In 
the competition between the AH and SH states at weak 
interaction and in finite uniform axial field we found that 
at half filling the SH state wins. The TRS breaking AH 
order paramater, however, appears at finite doping. We 
hope that the realization of time-reversal symmetric ax- 
ial field and its tunability over a wide range in real^I and 
artificial graphene^ 2 - will make these effects soon visible 
in experiments. 
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In the paper we have presented numerical evidence 
(for spinless fermions) that in the presence of axial mag- 
netic fields the TRS can be spontaneously broken at weak 
second-nearest-neighbor repulsion, so that the system en- 
ters into the AH insulator phase. Below we provide some 
additional numerical results in support of our claim. 
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strength of interaction for the insulation in the absence 
of axial fields. 




FIG. 6: Variation of the TRS-breaking AH order with the 
next-nearest-neighbor interaction (V) at zero axial field. In- 
set: The variation of AH order in the entire system on A and 
B sub-lattices for V — 3, 2, 1.5, 1.3, reads from top to bottom. 



FIG. 7: Self-consistent solution of the AH order on A (left) 
and B (right) sublattices of a honeycomb lattice of 600 sites, in 
the presence of uniform axial magnetic field with b — 0.035&0 
(top), 0.03&o (bottom). The red, black, blue, magenta, green 
curves correspond to V = 1.5, 1.27, 1.0, 0.75, 0.5, respectively. 



A. Zero field criticality 

Let us first present the numerical solution of AH OP 
for zero axial field and estimate the critical strength of 
interaction (V c ) for insulation. In the paper we have 
mentioned that in the absence of axial fields, the next- 
nearest-neighbor component of the Coulomb interaction 
(V) in graphene needs to be sufficiently large to develop 
a TRS-broken state. A nonzero, and fairly uniform self 
consistent solution of the AH OP, computed on a 600 site 
honeycomb lattice can only be found for V > 1.27, as 
shown Fig. [6] The magnitude of the OP on both sublat- 
tices is equal. For V < 1.27 the OP vanishes everywhere 
in the system, which is an artifact of the self-consistent 
Fock approximation. The value of V c does not depend 
on system's size, when it contains more than 600 lat- 
tice points. We can then define V = 1.27 as the critical 



B. Pseudo magnetic catalysis 

In the paper, we have presented the distribution of 
the AH OP on A and B sub-lattices for a particular 
strength of the uniform axial magnetic field. However, we 
have performed the numerical analysis for various other 
strengths of the uniform axial magnetic field. In Fig. 
we present the AH OP on two sub-lattices for two differ- 
ent strengths of uniform axial magnetic field. 

We have also searched for the self-consistent solution of 
AH order in the presence of a bell-shaped axial magnetic 
field, localized around the center of the system. In the 
paper we have presented the self-consistent solution of 
the AH OP on the sublattice A, for a particular choice of 
the total axial flux passing through the system. However, 
the same analysis has been performed for other values 



6 



FIG. 8: Self consistent AH OP on the honeycomb lattice of 
726 sites, in the presence of a nonuniform magnetic field with 
total axial flux $ totai = 9.42$ (left), 10.99$ (right). The 
strength of the interaction reads as V = 1.27, 1, 0.75, 0.5 from 
top to bottom. 



Finite size effects in non-uniform condensation 
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of the total axial flux penetrating the system, and the 
results are shown in Fig. |SJ 
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FIG. 9: AH OP in the vicinity of the localized flux as a func- 
tion of the system size for V = 1.0, and total flux 7.85 < J>o. 
Here we have normalized the OP with respect to the maxi- 
mum one, obtained on a 726 site honeycomb lattice. Inset: 
variation of the OP far away from the localized flux. 



FIG. 10: Variation of the the real component of the AH 
OP in the vicinity of the localized flux with the system size, 
for V = 1.0, and when a total flux 7.85<l?o. Inset: the same 
quantity far from the localized flux. 

In the paper, we have shown that in the presence of 
a bell-shaped axial magnetic field, localized around the 
center of the system, the AH OP dominantly develops 
in the vicinity of the localized axial flux, when the inter- 
action is sub-critical. Far away from the localized field, 
the OP disappears. Fig. |H] shows that such behavior 
is generic, and the OP very much saturates in the bulk 
of the honeycomb lattice of 726 sites, whereas that near 
the boundary of the system gradually disappears as the 
system size is increased. Furthermore, we have noticed 
that the AH OP in the presence of a non-uniform axial 
magnetic field is accompanied by a real component, for 
sub-critical interactions. In Fig. [TU] shows that appear- 
ance of such real component appears to be a finite size 
effect. 



